PRECONDITIONED HSS METHOD FOR FINITE ELEMENT 
APPROXIMATIONS OF CONVECTION-DIFFUSION EQUATIONS* 

ALESSANDRO RUSSQt AND CRISTINA TABLING POSSIOt 

Abstract. A two-step preconditioned iterative method based on the Hemiitian/Skew-Hermitian 
splitting is applied to the solution of nonsymnietric linear systems arising from the Finite Element 
- - . approximation of convection-diffusion equations. The theoretical spectral analysis focuses on the 

OO ' case of matrix sequences related to FE approximations on uniform structured meshes, by referring 

^ ^ ' to spectral tools derived from Toeplitz theory. In such a setting, if the problem is coercive, and the 

C 3 ' diffusive and convective coefficients are regular enough, then the proposed preconditioned matrix 

^Nj , sequence shows a strong clustering at unity, i.e., a superlinear preconditioning sequence is obtained. 

Under the same assumptions, the optimality of the PHSS method is proved and some numerical 
experiments confirm the theoretical results. Tests on unstructured meshes are also presented, showing 
the some convergence behavior. 
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^^ , 1. Introduction. The paper deals with the numerical solution of linear systems 

r^ • arising from the Finite Element approximation of the elliptic convection-diffusion 

C^ I problem 



div f — a(x)VM -1- /3(x)uj = /, x e fi, (i t\ 

u\an = 0. 

> 

f~^ ' We apply the two-step iterative method based on the Herniitian/Skew-Hermitian split- 

0^ , ting (HSS) of the coefficient matrix proposed in [3] for the solution of nonsymmetric 

\l_ ' linear systems whose real part is coercive. The aim is to study the preconditioning 

effectiveness when the Preconditioned HSS (PHSS) method is applied, as proposed 



C^ ' in [5]. According to the PHSS convergence properties, the preconditioning strategy 



for the matrix sequence {A„(a,/3)} can be tuned only with respect to the stiffness 
^^ I matrices. This allows to adopt the same preconditioning proposal just analyzed in 

the case of Finite Difference (FD) and Finite Element (FE) approximations of the 
. ^ _ diffusion problem in [HllTnilinillllllllllll], and, more recently, in the case of FD 

^ . approximations of (jl.ip . 

?H ' More precisely, we consider the preconditioning matrix sequence {P„(a)} defined as 

^ ■ Pr,{a) = Di^^{a)A^{l,0)Di^^{a), where D^{a) = diag(A„(a,0))diag-i(^„(l,0)), i.e., 

the suitable scaled main diagonal of An{a, 0). In such a way, the solution of the linear 
system with coefficient matrix A„(a, f3) is reduced to computations involving diago- 
nals and the matrix A„(l, 0). In the case of uniform structured meshes the latter task 
can be efficiently performed by means of fast Poisson solvers, among which we can list 
those based on the cyclic reduction idea (see e.g. JSKTUIIIS]) and several specialized 
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2 A. RUSSO AND C. TABLING POSSIO 

multigrid methods (see e.g. [01[T5] V 

Our theoretical analysis focuses on the case of matrix sequences {v4„(a, f3)} related to 
FE approximations on uniform structured meshes, since the powerful spectral tools 
derived from Tocplitz theory [SI [71 [TSl [IS] greatly facilitate the required spectral anal- 
ysis. In such a setting, under proper assumptions on a(x) and /3(x), we prove the 
optimality of the PHSS method. In our terminology (see 2 ) this means that the 
PHSS iterations number for reaching the solution within a fixed accuracy can be 
bounded from above by a constant independent of the dimension n — n{h). 
Nevertheless, the numerical experiments have been performed also in the case of 
unstructured meshes, with negligible differences in the PHSS method performances. 
In such cases, by coupling the PHSS method with standard Preconditioned Conju- 
gate Gradient (PCG) and Preconditioned Generalized Minimal Residual (PGMRES) 
method, our proposal makes only use of matrix vector products (for sparse or even 
diagonal matrices) and of a solver for the related diffusion equation with constant 
coefficient. 

The underlying idea is that whenever P„(a) results to be an approximate factoriza- 
tion of yl„(a,/3), the computational bottleneck lies in the solution of linear systems 
with coefficient matrix v4„(l, 0). Thus, the main effort in devising efficient algorithms 
must be devoted to this simpler problem. Moreover, to some extent, the approxima- 
tion schemes in the discretization should take into account this key step to give rise 
to a linear algebra problem less difficult to cope with. 

The outline of the paper is as follows. In section [2| we report a brief descrip- 
tion of the FE approximation of the convection-diffusion equation, while in section [3| 
we summarize the definition and the convergence properties of the Hcrmitian/Skew- 
Hermitian (HSS) method and of its preconditioned formulation (PHSS method). Sec- 
tion [4| analyzes the spectral properties of the matrix sequences arising from FE ap- 
proximations of the considered convection-diffusion problem and reports the precondi- 
tioner definition. Section[5|is devoted to the theoretical analysis of the preconditioned 
matrix sequence spectral properties in the case of structured uniform meshes. In sec- 
tion [5] several numerical experiments illustrate the claimed convergence properties 
and their extension in the case of other structured and unstructured meshes. Lastly, 
section [7| deals with complexity issues and perspectives. 

2. Finite Element approximation. Problem (jl.ip can be stated in variational 
form as follows: 

{find u G Hi {Q) such that 
^ faVu • Vv? - /3 • V(/7 w") = /f^ /(/3 for all Lp e H^in) '-^■^'* 

where Hllfl) is the space of square integrable functions, with L^ weak derivatives 
vanishing on dfl. We assume that fi is a polygonal domain and we make the following 
hypotheses on the coefficients 

aeC^iU), with a(x) > ao > 0, 

/3 € C^m, with div/3 > pointwise in n, (2.2) 

The previous assumptions guarantee existence and uniqueness for problem (j2.1|) . 
For the sake of simplicity, we restrict ourselves to linear finite element approximation 
of problem (|2.ip . To this end, let 7^ = {K} be a usual finite element partition of fl 
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into triangles, with hfc — diani(iir) and h — max^ hx- Let Vh C HQ{fl) be the space 
of hnear finite elements, i.e. 

Vh — {(fih : il —^ M. s.t. (fill is continuous, </''i|k i^ linear, and ^h\dn = 0}- 

The finite element approximation of problem (j2.ip reads: 

I find Uh S Vh such that 



/a(« 



\/uh ■ ^^h ~ ■ Vv3h Uh) ^ /j-j ffh for all iph &Vh. 



(2.3) 



For each internal node i of the mesh 7^, let (pi G Vh be such that ipi{node i) = 1, 
and (/9i(node j) = if i ^ j. Then, the collection of all ipi's is a base for Vh- We 
will denote by n{h) the number of the internal nodes of Th, which corresponds to the 
dimension of Vh ■ Then, we write Uh as 

n(h) 

U/i = ^ Uiipt 
i=l 

and the variational equation (j2.3p becomes an algebraic linear system: 

V" I / aVipj ■ Vip, - Vip, ■ (3 ipj) Uj = / /(/3j, i = 1, . . . ,n(/i). (2.4) 

The aim of this paper is to study the effectiveness of the proposed Preconditioned 
HSS method applied to the quoted nonsymmetric linear systems (|2.4[) , both from the 
theoretical and numerical point of view. 

3. Preconditioned HSS method. In this section we briefly summarize the def- 
inition and the relevant properties of the Hermitian/Skew-Hermitian (HSS) method 
formulated in [3] and of its extension in the case of preconditioning as proposed in [5] . 
The HSS method can be applied whenever we are looking for the solution of a linear 
system v4„x = b where An G C"'*" is a non singular matrix with a positive definite 
real part and x, b belong to C" . Several applications in scientific computing lead to 
such kind of linear problems and, typically, the matrix An is also large and sparse, as 
in the case of FD or FE approximations of (|l.ip . 

More in detail, the HSS method refers to the unique Hermitian/Skew-Hermitian split- 
ting of the matrix An as 

A„=Re(A„)+iIm(A„), i^ = -1 (3.1) 

where 

Re(A„)^ ^"+^" and Im(A„) ^ ^" ~.^" 

are Hermitian matrices by definition. 

In the same spirit of the ADI method [11] , the quoted splitting allows to define the 

Hermitian/Skew-Hermitian (HSS) method [3], as follows 



(a/ + Re(yl„))x'=+^ = (a/ - ilm(yl„))x''' + b 
(a/ + iIm(A„))x'=+i = (a/-Re(yl„))x'^+^ +b 



k+l - (r^T _Tir.( A \\^k+\ , 1- (3-2) 
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with a positive parameter and x° given initial guess. 

Beside the quoted formulation as a two-step iteration method, the HSS method can 

be reinterpreted as a stationary iterative method whose iteration matrix is given by 

M{a) = [al + iIm(A„))"^ {al - Rc(A„)) [al + Rc(A„))"^ {al - iIm(A„)) (3.3) 

and whose convergence properties are only related to the spectral radius of the Her- 

mitian matrix [al — Re(^„)) {al + Re(A„)) , which is unconditionally bounded by 

1 provided the positivity of a and of Re(^„) [3 . 

Indeed, the rate of convergence can be unsatisfactory for large values of n in the case 

of PDEs applications, as for instance (jl.ip . so that the preconditioned formulation of 

the method proposed in [S] is more profitable. 

Let Pn be a Hermitian positive definite matrix. The Preconditioned HSS (PHSS) 

method can be defined as 



(a/ + P-iRe(A„))x'=+^ = (a/-P-iiIm(yl„))x''^ + p-ib 
{al + p-iiIm(A„)) x'-^+i = [al - p-iRe(A„)) x'=+^ + p-ib 



(3.4) 



Notice that the proposed method differs from the HSS method applied to the matrix 
P^^An since P,7^Re(A„) and P^^Im(A„) are not the Hcrmitian/Skew-Hermitian 
sphtting of P^^An. 

Let \{X) denote the set of the eigenvalues of a square matrix X; the convergence 
properties of the PHSS method exactly mimic those of the previous HSS method, as 
claimed in the theorem below. 

Theorem 3.1. 5J Let An G C"^" he a matrix with positive definite real part, a 
be a positive parameter and let Pn G C"^" be a Hermitian positive definite matrix. 
Then the iteration matrix of the PHSS method is given by 

M{a) = {al + ip-ilm(A„))"' {al - p-iRe(A„)) {al + p-iRe(A„))^' 
(a/-ip-ilm(A„)), 

its spectral radius g{AI{a)) is bounded by 

a - Xi 



<7{a) = max 

AieA(P^^Ro(A„)) 



a + Xi 



< 1 for any a > 0, 



i.e., the PHSS iteration is unconditionally convergent to the unique solution of the 
system A„x = b. Moreover, denoting by k = Ai„ax(P,r^R,'3(v4„))/Aniin(P,r^R'6(A„)) 
the spectral condition number (namely the Euclidean (spectral) condition number of 
the symmetrized matrix), the optimal a value that minimizes the quantity cr(a) equals 

a* = A/A„,in(P„-iRc(A„))A„,ax(Pn 'Re(A„)) and a{a*)^^^^^^. 

Thus, the unconditional convergence property holds also in the preconditioned for- 
mulation of the method and the convergence properties are related to the spectral 
radius of 

[al p-^/'MAn)P-'/') [al + p-i/2Re(A„)p-i/2 

where the optimal parameter a is the square root of the product of the extreme 
eigenvalues of P,7^Re(yl„). 
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It is worth stressing that the PHSS method in (|3.4p can also be interpreted as the 
original iteration (|3.2p where the identity matrix is replaced by the preconditioner P„, 



(aP„ + Re(A„))x'=+5 = (aP„ -ilm(^„))x'= + b 

(aP„ + iIm(AO)x'=+i = (aP„-Re(A0)x'=+5+b ^'' 

Clearly, the last formulation is the most interesting from a practical point of view, 
since it does not involve any inverse matrix and it easily allows to define an inexact 
formulation of the quoted method: in principle, at each iteration, the PHSS method 
requires the exact solutions with respect to the large matrices aP„ + Re(A„) and 
aPn + ilm(yl„). This requirement is impossible to achieve in practice and an inexact 
outer iteration is computed by applying a Preconditioned Conjugate Gradient method 
(PCG) and a Preconditioned Generalized Minimal Residual method (PGMRES), with 
preconditioner P„, to the coefficient matrices aP„ + Re(yl„) and aPn + i Im(A„), 
respectively. Hereafter, we denote by IPHSS method the described inexact PHSS 
iterations. 

The accuracy for the stopping criterion of these additional inner iterative procedures 
must be chosen by taking into account the accuracy obtained by the current step of 
the outer iteration (see [S] H] and section [S] for some remarks about this topic) . The 
most remarkable fact is that the PHSS and IPHSS methods show the same conver- 
gence properties, though the computational cost of the latter is substantially reduced 
with respect to the former. 

In the IPHSS perspective the spectral properties induced by the preconditioner P„ 
are relevant not only with respect to the IPHSS rate of convergence, i.e., the outer 
iteration, but also with respect to the PCG and PGMRES ones. So, the spectral 
analysis of the matrix sequence {P,7^Im(yl„)} becomes relevant exactly as the spec- 
tral analysis of the matrix sequence {P^^Ke{An)}. 

Lastly, it is worth stressing that a deeper insight of the HSS/PHSS convergence prop- 
erties with respect to the skew-Hermitian part pertains to the framework of multi- 
iterative methods |14] , among which multigrid methods represent a classical example. 
Typically, a multi-iterative method is composed by two, or more, different iterative 
techniques, where each one is cheap and potentially slow convergent. Nevertheless, 
these iterations have a complementary spectral behavior, so that their composition 
becomes fast convergent. In fact, the PHSS matrix iteration can be reinterpreted as 
the composition of two distinct iteration matrices and the strong complementarity 
of these two components makes the contraction factor of the whole procedure much 
smaller than the contraction factors of the two distinct components. In particular, the 
skew-Hermitian contributions in the iteration matrix can have a role in accelerating 
the convergence. The larger is /3(x), i.e., the problem is convection dominated, the 
more the FE matrix A„ departs from normality. Thus, a stronger "mixing up effect" 
is observed and the real convergence behavior of the method is much faster compared 
with the forecasts of Theorem 13. II See [S] for further details. 



4. Preconditioning strategy. According to the notations and definitions in 
Section [21 the algebraic equations in ()2.4p can be rewritten in matrix form as the 
linear system 

An{a,(3)x = h, 
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with 

A„(a,/3) = e„(a) + *„(/3)eM"^", n^n{h), (4.1) 

where 6n(a) and ^n(/3) represent the approximation of the diffusive term and ap- 
proximation of the convective term, respectively. More precisely, we have 

{Qn{a))^,J = / aWip^WlpJ (4.2) 

(vI/„(/3)),,, =- /'(V^,-/3)<^„ (4.3) 

where suitable quadrature formula are considered in the case of non constant coeffi- 
cient functions a and (3. 

Thus, according to (|3.ip . the Hcrmitian/skew-Hermitian decomposition of A„(a, (3) is 
given by 

Re(A„(a,/3)) = e„(a) + Re(*„(/3)), (4.4) 

iIm(A„(a,/3)) = ilm(*„ (/?)), (4.5) 



where 



Re(vI/„(/3)) = i(*„(/3) 4- *^(/3)) = i?„(;9). 



iIm(*„(/3)) = *„(/3)-£;„(/3), 

since by definition, the diffusion therm 8n(a) is a Hermitian matrix and does not 

contribute to the skew-Hermitian part of A„(a,/3). Notice also that £'„(/3) = if 

div(/3) = 0. 

Clearly, the quoted HSS decomposition can be performed on any single elementary 

matrix related to T/j by considering the standard assembling procedure. 

By construction, the matrix Re(A„(a, /3)) is symmetric and positive definite whenever 

Amin(0n(a)) > p{En{(3)) ■ ludccd, without the condition div(/3) > 0, the matrix £'„(/?) 

does not have a definite sign. 

Moreover, the Lemma below allows to obtain further information regarding such a 

structural assumption. 

Lemma 4.1. Let {i?„(/3)} be the matrix sequence defined as 

i?«(/9) = ^(*n(/3)-fC(/3))- 



Under the assumptions in (j2.2p . then it holds 

\\En0)h < l|-E„(/3)||oo < Ch^, 

with C absolute positive constant only depending on P{:x.) and fi. 

Proof. By applying the Green formula, it holds that for any i, j = 1, . . . ,n{h) 

div(/3) • ipi^pj 



2./0 



2 XI / div(;3)(^,v?j 



KCS{ipi)nS{^j) 
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with respect to the related mesh 7/j — {K\ and where S{ipk) denotes the support of 
basis element ip^. on 7/j. Thus, we have 



qh^ 



\En0)U.j < -sup div(^) Y^ / \ip,ipj\ < -sup div(/3) 

since \ipk\ < 1 for any k and 

where h is the finesse parameter of the mesh 7^ and q equals the maximum number 
of involved mesh elements with respect to the mesh sequence {Th}. 
Lastly, since En{P) is a Hermitian matrix, it holds 



\En0)h < \\En0)\\oo < D . . max |£;„(/3)|.j < ^Dsup |div(/3')|q/i' 



2 



;j = l,...,n(h) 4 

where D denotes the maximum number of nonzero entries on the rows of £"„ (/?) . D 



Remark 4.2. The claim of Lemma 14.11 holds whenever a quadrature formula 
with error 0{h^) is consider for approximating the integrals involved in (|4.3p . 



Moreover, in the special case of a structured uniform mesh on il = (0, 1)^ as con- 
sidered in the next section, under the assumptions of Lemma |4. II and a(x) > oq > 0, 
we have 6„(a) > c/i^/„ (with c absolute positive constant), so that Re(^„(a, /3)) > 
(c— C)h'^In- Here, we are referring to the standard ordering relation between Hermi- 
tian matrices, i.e., the notation X >Y , with X and Y Hermitian matrices, means that 
X — Y is nonnegative definite. Thus, under the assumption that |div(/3)| is smaller 
than a positive suitable constant, it holds that Re(A„(a,/?)) is real, symmetric, and 
positive definite. 

However, the main drawback is due to ill-conditioning, since the condition number is 
asymptotic to /i^^, so that preconditioning is highly recommended. 
By referring to a preconditioning strategy previously analyzed in the case of FD ap- 
proximations XL JJli -2Qj 22lj 23, or FE approximations [21^ of the diffusion equation 
and recently applied to FD approximations [5] of (jl.ip . we consider the precondition- 
ing matrix sequence {Pn{a)} defined as 

Pn{a) = Dl{a)Ar,{l,Q)Dl{a) (4.6) 

where Dn{a) — diag(yl„(a, 0))diag^^(A„(l, 0)), i.e., the suitable scaled main diagonal 
of A„(a,0) and clearly yl„(a,0) equals 6„(a). 

In such a way, the solution of the linear system A„(a, f3)\i = f is reduced to computa- 
tions involving diagonals and the matrix A„(1,0). In the case of uniform structured 
mesh this task can be efficiently performed by considering fast Poisson solvers, such 
as those based on the cyclic reduction idea (see e.e. OIiniES]) and several specialized 
multigrid methods (see e.g. [T3lfT8]). 

It is worth stressing that the preconditioner is tuned only with respect to the diffu- 
sion matrix Qn{a) owing to the PHSS convergence properties highlighted in section 
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(s,t+1) (s+1,t+1) 




(s-1,t) / (s,t) / (s+1,t) 



(s-1,t-1) (s,t-1) 



Fig. 5.1. Uniform structured mesh giving rise to a Toeplitz matrix in the constant coefficient 
case a(x) = 1. 

[3l Indeed, the PHSS shows a convergence behavior mainly depending on the spec- 
tral properties of the matrix P~^(a)Re{An{a,P)). Nevertheless, the skew-Hermitian 
contribution may play a role in speeding up considerably the convergence (see [5] for 
further details). 

Hereafter, we denote by {A„(a,/?)}, n — n{h) the matrix sequence associated to a 
family of mesh {T^}, with decreasing finesse parameter h. As customary, the whole 
preconditioning analysis will refer to a matrix sequence instead to a single matrix, 
since the goal is to quantify the difficult of the linear system resolution in relation to 
the accuracy of the chosen approximation scheme. 

5. Spectral analysis and clustering properties in the case of structured 
uniform meshes. In the present section we analyze the spectral properties of the 
preconditioned matrix sequences 

{P~\a)ReiA„{aJ))} and {p-i(a)Im(A„(a, /3))} 

in the special case of J7 = (0, 1)^ with a structured uniform mesh as in Figure [??T1 by 
using the spectral tools derived from Toeplitz theory [SI [71 [TSl [Hj . The aim is to prove 
the optimality of the PHSS method, i.e., the PHSS iterations number for reaching the 
solution within a fixed accuracy can be bounded from above by a constant independent 
of the dimension n = n{h). A more in depth analysis is also considered for foreseeing 
the IPHSS method convergence behavior. 
We make reference to the following definition. 

Definition 5.1. |26| Let {An} be a sequence of matrices of increasing dimensions 
n and let g he a measurable function defined over a set K of finite and positive Lehesgue 
measure. The sequence {^n} is distributed as the measurable function g in the sense 
of the eigenvalues, i.e., {A„} ~;^ g if, for every F continuous, real valued and with 
bounded support, we have 

hm i^F(A,(A„)) = -^ I F{e{s))ds, 

where Aj(A„), j = 1, . . . ,n, denote the eigenvalues of An. 

The sequence {A„} is clustered at p if it is distributed as the constant function 

g{x) = p, i.e., for any £ > 0, 7^ {i | Ai(j4„) (f. (p — £,p + e)} = o{n). The sequence 
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{A„} is properly (or strongly) clustered at p if for any e > the number of the eigen- 
values of An not belonging to {p^e,p+e) can be bounded by a pure constant eventually 
depending on e, but not on n. 

First, we analyze the spectral properties of the matrix sequence {f ~^(a)Re(A„(a, (3))} 
that directly influences the convergence behavior of the PHSS method. We refer to a 
preliminary result in the case in which the convection term is not present. 

Theorem 5.2. [51] Let {A„(a,0)} and {Pn{a)} be the Hermitian positive defi- 
nite matrix sequences defined according to (j4.ip and (|4.6p . If the coefficient a(x) is 
strictly positive and belongs to C^(f2), then the sequence {Pn^{a)An{a, 0)} is properly 
clustered at 1. Moreover, for any n all the eigenvalues of P~^{a)An{a,0) belong to 
an interval [d, D] well separated from zero [Spectral equivalence property] . 

The extension of this claim in the case of the matrix sequence {Re(A„(a,/3))} with 
Re(v4„(a, /?)) y^ On{ci) can be proved under the additional assumptions of Lemma 14. II 

Theorem 5.3. Let {Re(A„(a, /?))} and {Pn{a)} be the Hermitian positive definite 
matrix sequences defined according to (|4.4p and (j4.6p . Under the assumptions in (j2.2p . 
then the sequence {P^^(a)Re(A„(a. /3))} is properly clustered at 1. Moreover, for any 
n all the eigenvalues o/P^^(a)Re(A„(a,/3)) belong to an interval [d, D] well separated 
from zero [Spectral equivalence property] . 

Proof. The proof technique refers to a previously analyzed FD case [22] and it is 
extended for dealing with the additional contribution given by En{(3). 
First, we consider the spectral equivalence property. Due to a similarity argument, we 

analyze the sequence {e-i(l)(e;(a) + -E;(/3))}, where X* = DfJ{a)XDnHa) and 
D„(a) = diag(9„(a))diag~"^(0„(l)). According to the assumptions and Theorem 7 
in [21j , we have the following asymptotic expansion 



so that 



e;(a) = e„(l) + h^Frnia) + oih-')Gnia) (5.1) 



e-i(l)e;(a) = /„ + e-\l){h^Fn{a) + o{h^)Gn{a)). 



Taking into account the order of the zeros of the generating function of the Toeplitz 
matrix 8„(1), we infer that there exists a constant ci so that ||6^^(1)||2 < Ci/i~^ 
[6l[T5]. Moreover, it also holds that 

\\E:m\2<-. — . (5.2) 

mmo a 
Therefore, by standard linear algebra, we have 

A„.ax(e-i(i)(e:(a) + K(^)) < \\e-Hi)mx^) + Ei0))\\2 

< \\Inh + h^G-\l)h\\Fn{a) + o(l)G„(a)||2 

+||e-i(i)||2||i?:(^)ll2 

< 1 + ci (\\Fn{a)\\2 + o(l)||G„(a)||2 + ^ 



mmn a 
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where ||i^n(a)||2 and ||G„(a)||2 are uniformly bounded since Fn{a) and Gn{a) are 
bounded symmetric sparse matrices by virtue of Theorem 7 in |21j . 
Conversely, a bound from below for Ainin(0,r^(l)(9*i(a) + E^{P)) requires a bit dif- 
ferent technique with respect to Theorem l5.21 owing to the presence of the convective 
term. Again by a similarity argument and by referring to the Courant-Fishcr Theo- 
rem, we have 

A„,:„(9..(i)-(9:» + Km - mi. -"Q"'""^-.'-); ^■■(/?))ft-'(')- 

>A,nin(e,-i(i)e:(a)) 

. x^e;^(l)i;*(;9)9;^(l)x 

-I- mm T^ 

x/o x-' X 

where 

\ mmn a 

as proved in [5T1 (25] and 

^^.^^ x^9,T^(l)g*(/3)9.T^(l)x ^ K,^{E10)) ^ Ch^ _ C 



x^o x^x Ainin(9„(l)) C2h'^ uoRn 0, C2minoa' 

being A,„iii(e„(l)) <C2h?. 

The proof of the presence of a proper cluster again makes use of a double similarity 
argument and of the asymptotic expansion in (|5.ip . so that we analyze the spectrum 
of the matrices 

Xn^In + e;^(l)(/i2F„(a) -f o{h^)Gn{a) + E10))Q~'^ {I) 

similar to the matrices 8,7^(1) (6* (a) -I- E^0)). 

As in the case of Theorem 15.21 we refer to the matrix U e M"^p, whose columns 
are made up by considering the orthonormal eigenvectors of 8n(l) corresponding to 
the eigenvalues Ai(6„(l)) > \e~^^h? , since we know [35] that for any sequence {£:„} 
decreasing to zero (as slowly as wanted) it holds 

#X,„ = 0{\e-'}), X,„ = {*|A.(e„(l)) < \e-'}h^}. 

Thus, we consider the following split projection 

U'^XnU ^Ip + Yp + Zp + Wp 

with Yp == h^DpU'^FnUDp, Zp = o{h^)DpU'^GnUDp and Wp = h'^DpU'^ El{l3)UDp, 
where Dp = diag ( A^ ^ ) G M^^^. The matrices Yp + Zp + Wp are of infinitesimal 

spectral norm since all the terms 1^, Zp 21] and Wp are too. In fact, by virtue of 
(|5.2p and the Dp matrix definition we have 

1 Ch"^ C 

\\Wph < \\Dp\\l\\E:m\2 < T-rrn— — ^ -^— M- 

\en \ h^ mmo a mmo a 
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Therefore, by applying the Cauchy interlacing theorem [T^], it directly follows that 
for any e > there exists n such that for any n > h (with respect to the the ordering 
induced by the considered mesh family) at least n — 0([£^^]) eigenvalues of the 
preconditioned matrix belong to the open interval (1 — e, 1 + e). □ 



Remark 5.4. The claim of Theorem 15.31 holds both in the case in which the 
matrix elements in (|4.2p and ()4.3p are evaluated exactly and whenever a quadrature 
formula with error 0{h^) is consider to approximate the involved integrals. 



The previous results prove the optimality both of the PHSS method and of the PCG, 
when applied in the IPHSS method for the inner iterations. However, still in the 
case of the IPHSS method, suitable spectral properties of the preconditioned matrix 
sequence {P^^(a)Im(A„(a, /3))} have to be proven with respect to the PGMRES ap- 
plication. 

Hereafter, even if we make direct reference to the spectral Toeplitz theory, we prefer 
to preliminary analyze the matrices by considering the standard FE assembling proce- 
dure. Indeed, the FE elementary matrices suggest the local domain analysis approach 
in a more natural way than in the FD case. More precisely, this purely linear algebra 
technique consists in an additive decomposition of the matrix in simpler matrices al- 
lowing a proper majorization for any single term (see also [5j|4]). 

Theorem 5.5. Let {Im(yl„(a,/3))} and {P„(a)} he the Hermitian matrix se- 
quences defined according to (|4.5p and ()4.6p . Under the assumptions in (|2.2|) . then 
the sequence {P^^(a)Im(A„(a, /3))} is spectrally bounded and properly clustered at 
with respect to the eigenvalues. 

Proof. Clearly, by referring to the standard assembling procedure, the Hermitian 
matrix Im(A„(a,/3)) can be represented as 

Im(A„(a,/3))= Y. Im(4f(a,/3)) 
KeTh 



with the elementary matrix corresponding to Im(j4„ (a,/3)) given by 



Im(A«(a,/3)) = -; 






-(712 -721) 


-(713-731) 


-721 





-(723-732) 


-731 


723 - 732 






Here 7^ — j{ipi,ipj) — j^iy^i ■ (3) ^j and the indices i,j refer to a local counter- 
clockwise numbering of the nodes of the generic element K E T^. 
With respect the desired analysis this matrix can be split, in turn, as 



Im(A^'(a,/3)) = -i((7: 



12-721)^^^^ 



(7 



23 



732 )Ag,^^ + (713 - 731 



)A^') , 



where 



A 



K.l 



0-10 
1 




A 



K.l 




0-1 
1 



A 



K,3 



0-1 



1 

,K,r 



By virtue of a permutation argument, the immersion of each matrix A^i'' , r = 1, . . . , 3, 
into the original matrix Im(A„(a, /?)) can be associated to a null matrix of dimension 
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,K.l ■ 



n{h) except for the 2 by 2 bfock given by A^^' in diagonal position. Therefore, we 
can refer to the same technique considered in 51]. More precisely, it holds that 



±i 



" 


-1 ■ 


< w 


[ 1 


u J 





l + /l2 -1 

-1 1 + K^ 



provided that w > {hyh? + 2) ^ and also for any constant 7 

< I7I w 



±7^ei =±71 



0-10 
1 




l + h^ -1 
-1 l + /i2 




\l\wR 



KA 



Thus, by linearity and positivity, we infer that 
'"- -^ - - (1 



±lva{Afi{aJ)) < I (I712 -72i|<'' + I723 - 732|i?,'I'' + |7i3 -73i|<'' 



<wh\\P\\ 



since 



R 



K,l 



-1 


-1 

1 + h^ 


" 



pK,2 








. 





2{l + h^) -1 -1 

-1 2(l + /i2) _i 

-1 -1 2{l + h^) 





l + h^ -1 
-1 l + h"^ 



pK,3 


' 1 + h^ 








-1 






-1 





l + h^ 



and |7ij - jj^\ < 2/i||/3||oo for any i,j. 

Therefore, with respect to the matrices of dimension n{h), we obtain the following 

key majorization 

±lm{A„{aJ)) < wh\\(3\\oo{cen{l) + 2h''l„), 

so that by virtue of LPOs and spectral equivalence properties proved in Theorem l5.3l 
it holds that 

±Im(A„(a,^)) < whWPWoo ( ^^e„(a) + 2h^lJ\ < wch\\^\\oo (Pn{a) + ch^In) ■ 

\ mmn a J \ / 

In such a way the claim can be obtained according to the very same proof technique 
considered in [5j. In fact, by considering the symmetrized preconditioned matrices 
and by referring to the Courant-Fisher theorem, we find 

±P„"^(a)Im(A„(a,^))P,7^(a) < wc/i||/3|U(/n + ch^P-Ha)) 

and the spectral analysis must focus on the spectral properties of the Hermitian matrix 
sequences {wh{In + ch'^Pj^^{a))}. 

By referring to Toeplitz spectral theory, the spectral boundeness property is proved 
since we simply have 

Amax(w/i(/„ + ch^p-\a))) <whil+"f—\ ^s + js-^ 
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independent of h, if we choose w — h^^ and h = he^^ for any given e > 0. 
Moreover, thanks to the strictly positivity of a(x) and for the same choice of w and 
h, we infer that 



c 



h^ 



<e + - A;i(e„(i))). 

nimfi a e 
Since #{s | A,((e„(l))) < h^/e} ^ 0(e"^) the proof is over. D 



Remark 5.6. The claim of Theorem 15.51 holds also in the case in which the 
matrix elements in (|4.3p are evaluated by applying any quadrature formula with error 
0{h^). 

To sum up, with the choice a — 1 and under the regularity assumptions of The- 
orems 15.31 and 15. 5[ we have proven that the PHSS method is optimally convergent 
(linearly, but with a convergence independent of the matrix dimension n{h) due to 
the spectral equivalence property) . In addition, when considering the IPHSS method, 
the PCG converges superlinearly owing to the proper cluster at 1 of the matrix se- 
quence {P~^(a)Re(^„(a, /3))}, that clearly induces a proper cluster at 1 for the matrix 
sequence {/ -I- P~^(a)Re(A„(a,/3))}. 

In an analogous way, the spectral boundeness and the proper clustering at of 
the matrix sequence {f ~^(a)Im(A„(a,/3))} allow to claim that PGMRES in the in- 
ner IPHSS iteration converges superlinearly when applied to the coefficient matrices 
{/ + iF-i(a)Im(A„)}. 

6. Numerical tests. Before analyzing in detail the obtained numerical results, 
we wish to give a general description of the performed numerical experiments and of 
some implementation details. The case of unstructured meshes is discussed at the end 
of the section, while further remarks on the computational costs are reported in the 
next section. 

Hereafter, we have applied the PHSS method, with the preconditioning strategy de- 
scribed in SectionlH to FE approximations of the problem (|1.1|) . First, we consider the 
case of a uniformly positive function and (3 function vector which are regular enough 
as required by Lemma 14.11 and Theorems 15.31 and 15.51 The domain of integration Q. 
is the simplest one, i.e., J7 — (0, 1)^ and we assume Dirichlet boundary conditions. 
Whenever required, the involved integrals have been approximated by means of the 
middle point rule (the approximation by means of the trapezoidal rule gives rise to 
analogous results) . 

As just outhned in Section [3l the preconditioning strategy has been tested by consid- 
ering the PHSS formulation reported in (|3.5|) and by applying the PCG and PGM- 
RES methods for the Hermitian and skew-Hermitian inner iterations, respectively. 
Indeed, in principle each iteration of the PHSS method requires the exact solutions 
with large matrices as defined in (j3.4p . which can be impractical in actual imple- 
mentations. Thus, instead of inverting the matrices a/ + P^^ {ayY{,e{An{a, j3)) and 
al -h P^i(a)Im(A„(a,^)), the PCG and PGMRES are applied for the solution of 
system with coefficient matrices aPn{a) + Re(^„(a, /?)) and aPn{a) -\- Im(A„(a,/9)), 
respectively, and with Pn{a) as preconditioner. 
Preliminary, the numerical tests have been performed by setting the tolerances re- 
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quired by the inner iterative procedures at the same value of the outer procedure 
tolerance. This allows a realistic check of the PHSS convergence properties in the 
case of the exact formulation of the method. Moreover, a significant reduction of the 
computational costs can be obtained by considering the inexact formulation of the 
method, denoted in short as IPHSS. In the IPHSS implementation of (13. 5p . the inner 
iterations in are switched to the (k + l)-th outer step if 

W'^J.PCgIU < (. , fc \\rj,PGMREs\\2 < (. , fc .p , n 

II II ^' ^ ' I •) II II ^' ^ 'I 1 \^ J 

IFklh iFklh 

respectively, where k is the current outer iteration, rj € (0,1), and where rj is the 
residual at the j-th step of the present inner iteration O [5] . The reported results 
refer to the case S — 0.9, that typically gives the best performances. 
The quoted criterion is effective enough to show the behavior of the inner and outer 
iterations for IPHSS. It allows to conclude that the IPHSS and the PHSS methods 
have the same convergence features, but the cost per iteration of the former is sub- 
stantially reduced as evident from the lower number of total inner iterations. 
It is worth stressing that more sophisticate stopping criteria may save a significant 
amount of inner iterations with respect to (j6.ip . In particular, the approximation 
error of the FE scheme could be taken into account to drive this tuning [Tj . 
Finally, mention has to be made to the choice of the parameter a. Despite the 
PHSS method is unconditionally convergent for any a > 0, a suitable tuning, ac- 
cording to Theorem 13.11 can significantly reduce the number of outer iterations. 
Clearly, the choice a = 1 is evident whenever a cluster at 1 of the matrix sequence 
{P,7^(a)Re(A„(a, (3))} is expected. In the other cases, the target is to approximatively 
estimate the optimal a value 



a = VA„,in(Pn-'Re(A„))An,ax(Pn-'Re(A„)) 

that makes the spectral radius of the PHSS iteration matrix bounded by (^{a*) = 
^/k~ 1/^/k+I, with K = XrnaxiPn^'Re{An))/XminiPn^'R.e{An)) spcctral Condition 
number of P~^Jie{An), namely the Euclidean (spectral) condition number of the 
symmetrized matrix. 

All the reported numerical experiments are performed in Matlab, with zero initial 
guess for the outer iterative solvers and stopping criterion ||r-fc||2 < 10~^||ro||2. 
No comparison is explicitly made with the case of the HSS method, since the obtained 
results are fully comparable with those observed in the FD approximation case [SJ [S] ■ 
In Table 16.11 we report the number of PHSS outer iterations required to achieve the 
convergence for increasing values of the coefficient matrix size n = n(h) when consid- 
ering the FE approximation with the structured uniform mesh reported in Figure [5T] 
and with template function a{x,y) = ai{x,y) = exp{x + y), f3{x,y) = [x y]'^ satis- 
fying the required regularity assumptions. The averages per outer step for PCG and 
PGMRES iterations are also reported (the total is in brackets); the values refer to the 
case of inner iteration tolerances that equals the outer iteration tolerance tol — 10"''. 
The numerical experiments plainly confirm the previous theoretical analysis in Sec- 
tion [5] In particular, we observe that the outer convergence behavior does not depend 
on the coefficient matrix dimension n — n{h). The same holds true with respect to 
the PCG inner iterations, while some dependency on n is observed with respect to 
the PGMRES inner iterations. More precisely, a higher number of PGMRES inner 
iterations is required in the first few steps of the PHSS outer iterations for increasing 
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Table 6.1 
Number of PHSS/IPHSS outer iterations and average per outer step for PCG and PGMRES 
inner iterations (total number of inner iterations in brackets). 





( 


i{x,y) = exp{x + y), /3( 


x,y) = X 


yV 




n 


PHSS 


PCG PGMRES 


IPHSS 


PCG 


PGMRES 


81 

361 

1521 

6241 

25281 


5 
5 
5 
5 
5 


1.6 (8) 2.4 (12) 
1.6 (8) 2.8 (14) 
1.6 (8) 3 (15) 
1.6 (8) 3.2 (16) 
1.6 (8) 3.6 (18) 


5 
5 
5 
5 
5 


1(5) 
1(5) 
1(5) 
1(5) 
1(5) 


1(5) 
1(5) 
2(10) 
2(10) 
2(10) 



Table 6.2 
Outliers analysis. 









a{x,y) = 


2xp{x + y), I3{x,y) = 


= [X yr 




n 


m,_ 


m-i- 


Pn\a)Rc 

Ptot 


(A„(a,/3)) 


'^max 


m_ 


P„-^(a)Im(A„(a,/3)) 

"1+ Ptot Ajnin 


-^max 


81 







3 


0% 
3% 


9.99e-01 


1.04e+00 




4 


0% 

4 9% 


-2.68e-02 


2.68e-02 


361 







4 


0% 
1% 


9.99e-01 


1.04e+00 




7 


0% 

7 3.8% 


-2.87e-02 


2.87e-02 


1521 







4 


0% 
0.26% 


9.99e-01 


1.044e+0 



9 


0% 
9 1.18% 


-2.93e-02 


2.93e-02 



n. Nevertheless, a variable tolerance stopping criterion for the inner iterations as 
devised in the IPHSS method is able to stabilize, or at least to strongly reduce, this 
sensitivity. 

The numerical results in Table 16.21 give evidence of the strong clustering properties 
when the previously defined preconditioner Pn(a) is applied. More precisely, for in- 
creasing values of the coefficient matrix dimension n, we report the number of outliers 
of P^^(a)Re(A„(a,/3)) with respect to a cluster at 1 with radius 6 = 0.1 (or S = 0.01): 
m_ is the number of outliers less then 1 — (5, m^ is the number of outliers greater 
then l + S, Ptot is the related total percentage. In addition, we report the minimal and 
maximal eigenvalue of the preconditioned matrices. The same information is reported 
for the matrices iP~^(a)Im(A„(a,/3)), but with respect to a cluster at 0. 

Despite the lack of the corresponding theoretical results, we want to test the 
PHSS convergence behavior also in the case in which the regularity assumption on 
a{x,y) in Theorems 15.31 and 15.51 are not satisfied. The analysis is motivated by 
favorable known numerical results in the case of FD approximations (see, for in- 
stance, [33J [231 ISj) or FE approximation with only the diffusion term [3Tj. More 
precisely, we consider as template the C^ function a{x,y) = a2{x,y) = e^+l^"^"' , 
the C° function a(x,y) = a^(x,y) = e^+l^~^/^l, and the piecewise constant function 
a{x,y) = 04(2;, y) = 1 if y < 1/2, 10 otherwise. 

The number of required PHSS outer iterations is listed in Table [6731 together with the 
averages per outer step for PCG and PGMRES inner iterations (the total is in brack- 
ets). Notice that in the case of the C^ or C" function the outer iteration number does 
not depend on the coefficient matrix dimension n = n{h). The same seems to be true 
with respect to the PCG inner iterations, while the PGMRES inner iterations show 
some influence on n. More precisely, this influence lies in a higher number of PGM- 
RES inner iterations in the first few steps of the PHSS outer iterations. Moreover, 
the considered variable tolerance stopping criterion for the inner iterations devised in 
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Table 6.3 
Number of PHSS/IPHSS outer iterations and average per outer step for PCG and PGMRES 
inner iterations (total number of inner iterations in brackets). 







a2{x 


,y), 0{x,y) -- 


= xyV 






n 


PHSS 


PCG 


PGMRES 


IPHSS 


PCG 


PGMRES 


81 


6 


2.2 (13) 


2.8 (17) 


6 


1(6) 


1(6) 


361 


6 


2.2 (13) 


3.2 (19) 


6 


1(6) 


2 (12) 


1521 


6 


2.2 (13) 


3.5 (21) 


6 


1(6) 


2 (12) 


6241 


6 


2.2 (13) 


4(24) 


6 


1(6) 


2(12) 


25281 


6 


2.2 (13) 


4.2 (25) 


6 


1(6) 


3(18) 






a-iix 


:»/): P{x,y) -- 


= xyV 






n 


PHSS 


PCG 


PGMRES 


IPHSS 


PCG 


PGMRES 


81 


7 


1.9 (13) 


2.6 (18) 


7 


1(7) 


1(7) 


361 


7 


2.2 (15) 


3(21) 


7 


1(7) 


1.7 (12) 


1521 


7 


2.2 (15) 


3.5 (24) 


7 


1.1 (8) 


2 (14) 


6241 


7 


2.3 (16) 


3.6 (25) 


7 


1.1 (8) 


2(14) 


25281 


7 


2.3 (16) 


4(28) 


7 


1.1 (8) 


2.1 (15) 



the IPHSS method is just able to reduce this sensitivity. 

These remarks are in perfect agreement with the outhers analysis of the matrices 
P~-^(a)Re(A„(a,/?)) and P~-^(a)Im(A„(a,/3)), with respect to a cluster at 1 and at 0, 
respectively, reported in Table WM with the same notations as before. 
Separate mention has to be made to the case of the piecewise continuous function 
a(x,y) = 04(2;, y). In fact, even if {P^^(a)Im(A„(a,/3))} could be supposed to be 
strongly clustered at 0, it is evident that {P'~^(a)Re(A„(o,/3))} is clustered at 1, but 
not in a strong way: the number of the outliers grows for increasing n, though their 
percentage is decreasing, in accordance with the notion of weak clustering (see Table 
I6.5p . Indeed, the number of outer iterations grows for increasing n as shown in Table 
16.61 and a deeper insight allows to observe that the major difficulty is, as expected, 
in the PCG inner iterations. The same behavior is observed also when varying the 
parameter a, in order to find the optimal setting (see Table [67f|) . 

Lastly, we want to test our proposal in the case of other structured and unstruc- 
tured meshes generated by triangle [24j with a progressive refinement procedure. The 
first meshes in the considered mesh sequences are reported in Figures 16.1116.31 
Tables 16.8116. lOl report the number of required PHSS/IPHSS iterations in the case 
of the previous template functions. Negligible differences in the PHSS/IPHSS outer 
iterations are observed for increasing dimensions n. Again, some dependency on n is 
observed with respect to the PGMRES inner iterations, due to an higher number of 
PGMRES inner iterations required in the first few steps of the PHSS outer iterations 
in relation to a more severe ill-conditioning. Clearly, a more sophisticated stopping 
criterion may probably reduce this sensitivity. 
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Table 6.4 
Outliers analysis. 









0-2(3:, y), 0{x,y) 


= [:. 


yV 




1 


n 


P„-^(a)Re(A„(a,/3)) 


p-^(a)Im(A„(a,/3)) | 


81 






1 

9 


1.2% 
11% 


9.97e-01 1.12e+00 




7 




7 


0% 
17% 


-4.32e-02 


4.32e-02 


361 






1 

11 


0.27% 
3% 


9.99e-01 1.12e+00 



15 




15 


0% 
8.3% 


-4.68e-02 


-4.68e-02 


1521 






1 

12 


6% 
0.79% 


9.99e-01 1.12e+00 




21 



21 


0% 
2.8% 


-4.78e-02 


4.78e-02 








«3(a;,2/), li{x,y) 


= [x 


yV 






n 


P„-^(a)Re(A„(a,/3)) 


p-^(a)Im(A„(a,/3)) | 


81 






1 

9 


1.2% 
11 % 


9.95e-01 1.16e+000 



6 




6 


0% 

14% 


-3.97e-02 


3.97e-02 


361 






1 
11 


0.28% 
3% 


9.97e-01 1.17e+00 



13 



13 


0% 
7% 


-4.31e-02 


4.31e-02 


1521 






1 

14 


0.07% 
0.92% 


9.98e-01 1.18e+00 




18 




18 


0% 
2.4% 


-4.40e-02 


4.40e-02 



Table 6.5 
Outliers analysis. 









a4,{x,y), 0{x,y) = 


[x y 


T 


1 


n 


p-\a)Ko{Ar,{aJ)) 


p-^(a)Im(A„(a,/3)) | 


81 


9 
9 


7 

9 


19% 
22% 


5.84e-01 2.09e+00 



1 




1 


0% 
2.5% 


-2.23e-02 2.23e-02 


361 


19 
19 


17 
20 


9.8% 
10.8% 


4.20e-01 2.97e+00 



3 



3 


0% 
1.6% 


-2.99e-02 2.99e-02 


1521 


39 
39 


37 
40 


5% 
5.2% 


2.78e-01 4.536+000 



6 



6 


0% 
0.79% 


-3.34e-02 3.34e-02 



Table 6.6 
Number of PHSS/IPHSS outer iterations and average per outer step for PCG and PGMRES 
inner iterations (total number of inner iterations in brackets). 







a4(x,y), P{x,y) 


= [x yf 






n 


PHSS 


PCG PGMRES 


IPHSS 


PCG 


PGMRES 


81 
361 
1521 


13 
20 
31 


3.5 (45) 2.2 (29) 
3.9 (78) 2.3 (47) 
4.3 (134) 2.4 (74) 


13 
20 
32 


1.9 (25) 

2.1 (41) 

2.2 (70) 


1(13) 
1(20) 
1.5 (47) 



Table 6.7 
Number of PHSS/IPHSS outer iterations and average per outer step for PCG and PGMRES 
inner iterations (total number of inner iterations in brackets) in the case of optimal a* values. 









a4,{x,y), I3{x,y) 


= [X yV 








n 


PHSS 


a* 


PCG PGMRES 


IPHSS 


a* 


PCG 


PGMRES 


81 
361 
1521 


12 
19 
29 


(1.068,1.12) 
(1.04,1.1656) 
(1.02,1.0736) 


3.8 (45) 2.2 (27) 
4.1 (77) 2.4 (45) 
4.5 (131) 2.4 (70) 


12 
18 
30 


(1.066,1.114) 
(1.088,1.1024) 
(1.0526, 1.16) 


2(24) 
2.1 (37) 
2.1 (64) 


1(12) 
1(18) 
1.4 (43) 
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Table 6.8 
Number of PHSS/IPHSS outer iterations and average per outer step for PCG and PGMRES 
inner iterations (total number of inner iterations in brackets) - meshes in Fig. \6.1l 







ai{os,y), P{x,y) = 


= [x yV 






n 


PHSS 


PCG PGMRES 


IPHSS 


PCG 


PGMRES 


41 

181 

761 

3121 

12641 

50881 


5 
5 
5 
5 
5 
5 


2.2 (11) 2.2 (11) 
2.2 (11) 2.6 (13) 
2.2 (11) 3 (15) 
2.2 (11) 3 (15) 
2.2 (11) 3.4 (17) 
2.2 (11) 3.8 (19) 


5 
5 
5 
5 
5 
5 


1(5) 
1(5) 
1(5) 
1(5) 
1(5) 
1(5) 


1(5) 
1(5) 
1.2 (6) 
2(10) 
2 (10) 
2(10) 






a,2{x,y), li{x,y) = 


= [x yV 






n 


PHSS 


PCG PGMRES 


IPHSS 


PCG 


PGMRES 


41 

181 

761 

3121 

12641 

50881 


6 
6 
6 
6 
6 
6 


2.2 (13) 2.7 (16) 
2.2 (13) 3 (18) 
2.2 (13) 3.3 (20) 
2.2 (13) 3.7 (22) 
2.2 (13) 4 (24) 
2.2 (13) 4.3 (26) 


6 
6 
6 
6 
6 
6 


1(6) 
1(6) 
1(6) 
1(6) 
1(6) 
1(6) 


1(6) 

1(6) 

2(12) 

2(12) 

2.2 (13) 

3(18) 






a-i{x,y), 0{x,y) = 


= [x yV 






n 


PHSS 


PCG PGMRES 


IPHSS 


PCG 


PGMRES 


41 

181 

761 

3121 

12641 

50881 


7 
7 
7 
7 
7 
7 


2.2 (15) 2.5 (17) 
2.2 (15) 2.8 (20) 
2.2 (15) 3.1 (22) 
2.4 (17) 3.6 (25) 
2.4 (17) 3.8 (27) 
2.4 (17) 3.9 (31) 


7 
7 
7 
7 
7 
8 


1(7) 
1(7) 
1.1 (8) 
1.1 (8) 
1.1 (8) 
1.1 (9) 


1(7) 
1(7) 
2(14) 
2 (14) 
2 (14) 
2.2 (18) 



Table 6.9 
Number of PHSS/IPHSS outer iterations and average per outer step for PCG and PGMRES 
inner iterations (total number of inner iterations in brackets) - meshes in Fig. 16.21 







aiix,y), 0{x,y) = 


= [x yV 






n 


PHSS 


PCG PGMRES 


IPHSS 


PCG 


PGMRES 


25 

113 

481 

1985 

8065 

32513 


5 
5 
5 
5 
5 
5 


2.2 (11) 2.2(11) 
2.2 (11) 2.4 (12) 
2.2 (11) 2.8 (14) 
2.2 (11) 3 (15) 
2.2 (11) 3.4 (17) 
2.2 (11) 3.6 (18) 


5 
5 
5 
5 
5 
5 


1(5) 
1(5) 
1(5) 
1(5) 
1(5) 
1(5) 


1(5) 

1(5) 

1(5) 

2(10) 

2(10) 

2(10) 






a2{x,y), f3{x,y) = 


= [x yV 






n 


PHSS 


PCG PGMRES 


IPHSS 


PCG 


PGMRES 


25 

113 

481 

1985 

8065 

32513 


6 
6 
6 
6 
6 
6 


2.2 (13) 2.5 (15) 
2.2 (13) 3 (18) 
2.2 (13) 3.2 (19) 
2.2 (13) 3.7 (22) 
2.2 (13) 4 (24) 
2.2 (13) 4.3 (26) 


6 
6 
6 
6 
6 
6 


1(6) 
1(6) 
1(6) 
1(6) 
1(6) 
1(6) 


1(6) 
1(6) 
2(12) 
2(12) 

2 (12) 

3 (18) 






a3(x,y), P{x,y) = 


= [X yf 






n 


PHSS 


PCG PGMRES 


IPHSS 


PCG 


PGMRES 


25 

113 

481 

1985 

8065 

32513 


6 

7 
7 
7 
7 
8 


2.2 (13) 2.5 (15) 
2 (14) 2.7 (19) 
2.1 (15) 3 (21) 

2.3 (16) 3.4 (24) 

2.4 (17) 3.8 (27) 
2.1 (17) 3.8 (30) 


6 

7 
7 
7 
7 
8 


1(6) 
1(7) 
1.1 (8) 
1.1 (8) 
1.1 (8) 
1.1 (9) 


1(6) 
1(7) 
1.8 (13) 
2 (14) 
2(14) 
2.1 (17) 
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Table 6.10 
Number of PHSS/IPHSS outer iterations and average per outer step for PCG and PGMRES 
inner iterations (total number of inner iterations in brackets) - meshes in Fig. 16. gl 







ai{x,y), P{x,y) = 


= xyV 






n 


PHSS 


PCG PGMRES 


IPHSS 


PCG 


PGMRES 


55 

142 

725 

1538 

7510 

15690 


5 
5 
5 
5 
5 
5 


2.2 (11) 2.2 (11) 
2.2 (11) 2.6 (13) 
2.2 (11) 3 (15) 
2.2 (11) 3 (15) 
2.2 (11) 3 (17) 
2.2 (12) 3.4 (18) 


5 
5 
5 
5 
5 
6 


1(5) 
1(5) 
1(5) 
1.2 (6) 
1.2 (6) 
1.2 (7) 


1(5) 
1(5) 
1(5) 
1.2 (6) 
1.8 (9) 
2(12) 






a2{x,y), I3ix,y) = 


= [x yV 






n 


PHSS 


PCG PGMRES 


IPHSS 


PCG 


PGMRES 


55 

142 

725 

1538 

7510 

15690 


6 
6 
6 
6 

7 
7 


2.2 (13) 2.7 (16) 
2.2 (13) 3 (18) 
2.2 (13) 3.3 (20) 
2.2 (13) 3.5 (21) 
2 (14) 3.7 (26) 
2 (14) 3.7 (26) 


6 
6 
6 
6 

7 
7 


1(6) 
1(6) 
1(6) 
1.2 (7) 
1.1 (8) 
1.1 (8) 


1(6) 

1(6) 

2(12) 

2(12) 

2(14) 

2(14) 






a3{x,y), P(x,y) = 


= [X yf 






n 


PHSS 


PCG PGMRES 


IPHSS 


PCG 


PGMRES 


55 

142 

725 

1538 

7510 

15690 


7 
7 
7 
7 
7 
8 


1.8 (13) 2.4 (17) 
2.2 (15) 2.8 (20) 
2.6 (18) 3.1 (22) 
2.6 (18) 3.4 (24) 
2.6 (18) 3.7 (26) 
2.4 (19) 3.6 (29) 


7 
7 
7 
7 
7 
8 


1(7) 
1(7) 
1.1 (8) 
1.1 (8) 
1.1 (8) 
1.1 (9) 


1(7) 
1(7) 
2 (14) 
2 (14) 
2 (14) 
2(16) 






Fig. 6.1. Structured meshes. 
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Fig. 6.3. Unstructured meshes. 

7. Complexity issues and perspectives. Lastly, we report some remarks 
about the computational costs of the proposed iterative procedure by referring to 
the optimality definition below. 

Definition 7.1. \2\ Let {A,„x,„ = b,„} he a given sequence of linear systems of 
increasing dimensions. An iterative method is optimal if 

1. the arithmetic cost of each iteration is at most proportional to the complexity 
of a matrix vector product with matrix Am , 

2. the number of iterations for reaching the solution within a fixed accuracy can 
he hounded from ahove by a constant independent of m. 



In other words, the problem of solving a linear system with coefficient matrix Am is 
asymptotically of the same cost as the direct problem of multiplying Am by a vector. 
Since we are considering the preconditioning matrix sequence {Pn{a)\ defined as 
Pn{a) = Ay'(a)A„(l,0)Ay'(a), where D,,{a) = diag(A„(a,0))diag-i(A„(l,0)), the 
solution of the linear system in (j2.4p with matrix A„{a, (3) is reduced to computations 
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involving diagonals and the matrix A„(l, 0). 

As well known, whenever the domain is partitioned by considering a uniform struc- 
tured mesh this latter task can be efficiently performed by means of fast Poisson 
solvers, among which we can list those based on the cyclic reduction idea (see e.g. 
[8l[T0l[25]) and several speciahzed multigrid methods (see e.g. O [18]). Thus, in 
such a setting, and under the regularity assumptions (|2.2|) . the optimality of the 
PHSS method is theoretically proved: the PHSS iterations number for reaching the 
solution within a fixed accuracy can be bounded from above by a constant indepen- 
dent of the dimension n = n{h) and the arithmetic cost of each iteration is at most 
proportional to the complexity of a matrix vector product with matrix A„(a,/3). 

Finally, we want to stress that the PHSS numerical performances do not get worse 
in the case of unstructured meshes. In such cases, again, our proposal makes only use 
of matrix vector products (for sparse or even diagonal matrices) and of a solver for 
the related diffusion equation with constant coefficient. To this end, the main effort 
in devising efficient algorithms must be devoted only to this simpler problem. 
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